{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "5aecd873",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.spatial import cKDTree\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# 设置真实周期（单位：天）\n",
    "true_period = 0.394\n",
    "\n",
    "# ========== STEP 1: 读取两个文件 ==========\n",
    "file_gr = \"/mnt/7b21f1e1-eb25-4cd5-bdb5-06d7d82fa253/Temp/force_photmetry/result50cm/asteroid_photometry_results.csv\"\n",
    "file_uviz = \"/mnt/7b21f1e1-eb25-4cd5-bdb5-06d7d82fa253/Temp/force_photmetry/result/asteroid_photometry_results.csv\"\n",
    "\n",
    "df_gr = pd.read_csv(file_gr)\n",
    "df_uviz = pd.read_csv(file_uviz)\n",
    "\n",
    "# 合并数据\n",
    "df_all = pd.concat([df_gr, df_uviz], ignore_index=True)\n",
    "df_all = df_all[df_all['calmag'].notna()].copy()\n",
    "df_all['timestamp'] = pd.to_datetime(df_all['obs_time'])\n",
    "df_all['mjd'] = (df_all['timestamp'] - df_all['timestamp'].min()).dt.total_seconds() / 86400.0\n",
    "df_all['phase'] = (df_all['mjd'] / true_period) % 1\n",
    "\n",
    "# ========== STEP 2: 绘制光变曲线 ==========\n",
    "band_order = ['u', 'v', 'g', 'r', 'i', 'z']\n",
    "colors = {'u': 'purple', 'v': 'blue', 'g': 'limegreen', 'r': 'red', 'i': 'orange', 'z': 'brown'}\n",
    "\n",
    "plt.figure(figsize=(9, 6))\n",
    "for band in band_order:\n",
    "    band_df = df_all[df_all['band'] == band]\n",
    "    if len(band_df) > 0:\n",
    "        plt.errorbar(band_df['phase'], band_df['calmag'], yerr=band_df['magerr'], fmt='o',\n",
    "                     label=f'{band} band', color=colors.get(band, 'gray'), alpha=0.7)\n",
    "plt.xlabel(\"Phase (folded by 0.394 days)\")\n",
    "plt.ylabel(\"Magnitude\")\n",
    "plt.title(\"Phase-folded Light Curves (g, r, u, v, i, z)\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.grid(True)\n",
    "plt.savefig(\"all.png\")\n",
    "plt.close()\n",
    "\n",
    "# ========== STEP 3: 颜色拟合分析函数 ==========\n",
    "def double_sine_model(phase, A1, A2, phi1, phi2, offset):\n",
    "    return A1 * np.sin(2 * np.pi * (phase + phi1)) + A2 * np.sin(4 * np.pi * (phase + phi2)) + offset\n",
    "\n",
    "def analyze_color_sigma_clipped(df_all, band1, band2, color_label, max_sep=0.02, sigma_clip=3.0):\n",
    "    df1 = df_all[df_all['band'] == band1][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band1}'}).reset_index(drop=True)\n",
    "    df2 = df_all[df_all['band'] == band2][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band2}'}).reset_index(drop=True)\n",
    "\n",
    "    tree = cKDTree(df2[['mjd']])\n",
    "    dists, idxs = tree.query(df1[['mjd']], distance_upper_bound=max_sep)\n",
    "    valid = idxs < len(df2)\n",
    "\n",
    "    df1_valid = df1[valid].copy()\n",
    "    df2_matched = df2.iloc[idxs[valid]].reset_index(drop=True)\n",
    "    df1_valid[f'mag_{band2}'] = df2_matched[f'mag_{band2}']\n",
    "    df1_valid['color'] = df1_valid[f'mag_{band1}'] - df1_valid[f'mag_{band2}']\n",
    "    df1_valid['phase'] = (df1_valid['mjd'] / true_period) % 1\n",
    "\n",
    "    # 去除 NaN/inf 并做 3σ 剔除\n",
    "    df1_valid = df1_valid[np.isfinite(df1_valid['color'])]\n",
    "    mean_color = df1_valid['color'].mean()\n",
    "    std_color = df1_valid['color'].std()\n",
    "#     df1_valid = df1_valid[np.abs(df1_valid['color'] - mean_color) < sigma_clip * std_color]\n",
    "\n",
    "    # 拟合\n",
    "    phase = df1_valid['phase'].values\n",
    "    color = df1_valid['color'].values\n",
    "\n",
    "    if len(phase) < 5:\n",
    "        print(f\"[{color_label}] 有效点过少，跳过拟合\")\n",
    "        return\n",
    "\n",
    "    p0 = [0.05, 0.02, 0.0, 0.0, np.mean(color)]\n",
    "    try:\n",
    "        popt, _ = curve_fit(double_sine_model, phase, color, p0=p0)\n",
    "        phase_fit = np.linspace(0, 1, 300)\n",
    "        color_fit = double_sine_model(phase_fit, *popt)\n",
    "\n",
    "        plt.figure(figsize=(7, 4))\n",
    "        plt.scatter(phase, color, alpha=0.7, label='Observed (clipped)', color='darkcyan')\n",
    "        plt.plot(phase_fit, color_fit, 'k-', lw=2, label='Double sine fit')\n",
    "        plt.xlabel(\"Phase (folded by 0.394 days)\")\n",
    "        plt.ylabel(color_label)\n",
    "        plt.title(f\"{color_label} Color Variation with Phase (3σ Clipped)\")\n",
    "        plt.grid(True)\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "        plt.savefig(f\"{color_label} Color Variation with Phase.png\")\n",
    "        plt.close()\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"[{color_label}] 拟合失败: {e}\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bb6451d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "4d9ef5de",
   "metadata": {},
   "outputs": [],
   "source": [
    "color_combinations = [\n",
    "    ('u', 'v', 'u - v'),\n",
    "    ('u', 'g', 'u - g'),\n",
    "    ('u', 'r', 'u - r'),\n",
    "    ('u', 'i', 'u - i'),\n",
    "    ('u', 'z', 'u - z'),\n",
    "\n",
    "    ('v', 'g', 'v - g'),\n",
    "    ('v', 'r', 'v - r'),\n",
    "    ('v', 'i', 'v - i'),\n",
    "    ('v', 'z', 'v - z'),\n",
    "\n",
    "    ('g', 'r', 'g - r'),\n",
    "    ('g', 'i', 'g - i'),\n",
    "    ('g', 'z', 'g - z'),\n",
    "\n",
    "    ('r', 'i', 'r - i'),\n",
    "    ('r', 'z', 'r - z'),\n",
    "    ('i', 'z', 'i - z'),\n",
    "]\n",
    "for b1, b2, label in color_combinations:\n",
    "    analyze_color_sigma_clipped(df_all, b1, b2, label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "8ebc9190",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_color_color(df_all, band1a, band1b, band2a, band2b, label_x, label_y, max_sep=0.02):\n",
    "    # 计算两个颜色量（差值）\n",
    "    df1 = df_all[df_all['band'] == band1a][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band1a}'}).reset_index(drop=True)\n",
    "    df2 = df_all[df_all['band'] == band1b][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band1b}'}).reset_index(drop=True)\n",
    "    df3 = df_all[df_all['band'] == band2a][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band2a}'}).reset_index(drop=True)\n",
    "    df4 = df_all[df_all['band'] == band2b][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band2b}'}).reset_index(drop=True)\n",
    "\n",
    "    tree1 = cKDTree(df2[['mjd']])\n",
    "    _, idxs1 = tree1.query(df1[['mjd']], distance_upper_bound=max_sep)\n",
    "    valid1 = idxs1 < len(df2)\n",
    "    df_c1 = df1[valid1].copy()\n",
    "    df_c1[f'mag_{band1b}'] = df2.iloc[idxs1[valid1]].reset_index(drop=True)[f'mag_{band1b}']\n",
    "    df_c1['color_x'] = df_c1[f'mag_{band1a}'] - df_c1[f'mag_{band1b}']\n",
    "\n",
    "    tree2 = cKDTree(df4[['mjd']])\n",
    "    _, idxs2 = tree2.query(df3[['mjd']], distance_upper_bound=max_sep)\n",
    "    valid2 = idxs2 < len(df4)\n",
    "    df_c2 = df3[valid2].copy()\n",
    "    df_c2[f'mag_{band2b}'] = df4.iloc[idxs2[valid2]].reset_index(drop=True)[f'mag_{band2b}']\n",
    "    df_c2['color_y'] = df_c2[f'mag_{band2a}'] - df_c2[f'mag_{band2b}']\n",
    "\n",
    "    # 尝试匹配两组颜色时间\n",
    "    tree_match = cKDTree(df_c2[['mjd']])\n",
    "    _, idxs_match = tree_match.query(df_c1[['mjd']], distance_upper_bound=max_sep)\n",
    "    valid = idxs_match < len(df_c2)\n",
    "\n",
    "    x = df_c1[valid]['color_x'].values\n",
    "    y = df_c2.iloc[idxs_match[valid]]['color_y'].values\n",
    "\n",
    "    # 绘图\n",
    "    plt.figure(figsize=(6, 5))\n",
    "    plt.scatter(x, y, alpha=0.7, color='slateblue')\n",
    "    plt.xlabel(label_x)\n",
    "    plt.ylabel(label_y)\n",
    "    plt.title(f\"{label_y} vs {label_x} (Color–Color)\")\n",
    "    plt.grid(True)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f\"{label_y} vs {label_x} (Color–Color).png\")\n",
    "    plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "daac7a08",
   "metadata": {},
   "outputs": [],
   "source": [
    "# g - r vs r - i\n",
    "plot_color_color(df_all, 'g', 'r', 'r', 'i', 'g - r', 'r - i')\n",
    "\n",
    "# u - g vs g - r\n",
    "plot_color_color(df_all, 'u', 'g', 'g', 'r', 'u - g', 'g - r')\n",
    "\n",
    "# v - z vs g - r\n",
    "plot_color_color(df_all, 'v', 'z', 'g', 'r', 'v - z', 'g - r')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b2924f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 增加拟合点数并加大图像分辨率以提高显示效果\n",
    "\n",
    "def analyze_color_sigma_clipped_v2(df_all, band1, band2, color_label, max_sep=0.02, sigma_clip=4.0):\n",
    "    df1 = df_all[df_all['band'] == band1][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band1}'}).reset_index(drop=True)\n",
    "    df2 = df_all[df_all['band'] == band2][['mjd', 'calmag']].rename(columns={'calmag': f'mag_{band2}'}).reset_index(drop=True)\n",
    "\n",
    "    tree = cKDTree(df2[['mjd']])\n",
    "    dists, idxs = tree.query(df1[['mjd']], distance_upper_bound=max_sep)\n",
    "    valid = idxs < len(df2)\n",
    "\n",
    "    df1_valid = df1[valid].copy()\n",
    "    df2_matched = df2.iloc[idxs[valid]].reset_index(drop=True)\n",
    "    df1_valid[f'mag_{band2}'] = df2_matched[f'mag_{band2}']\n",
    "    df1_valid['color'] = df1_valid[f'mag_{band1}'] - df1_valid[f'mag_{band2}']\n",
    "    df1_valid['phase'] = (df1_valid['mjd'] / true_period) % 1\n",
    "\n",
    "    # 去除 NaN/inf 并做 3σ 剔除\n",
    "    df1_valid = df1_valid[np.isfinite(df1_valid['color'])]\n",
    "    mean_color = df1_valid['color'].mean()\n",
    "    std_color = df1_valid['color'].std()\n",
    "    df1_valid = df1_valid [np.abs(df1_valid['color'] - mean_color) < sigma_clip * std_color]\n",
    "\n",
    "    # 拟合\n",
    "    phase = df1_valid['phase'].values\n",
    "    color = df1_valid['color'].values\n",
    "\n",
    "    if len(phase) < 5:\n",
    "        print(f\"[{color_label}] 有效点过少，跳过拟合\")\n",
    "        return\n",
    "\n",
    "    p0 = [0.1, 0.05, 0.0, 0.0, np.mean(color)]\n",
    "    try:\n",
    "        popt, _ = curve_fit(double_sine_model, phase, color, p0=p0)\n",
    "        phase_fit = np.linspace(0, 1, 500)  # 增加拟合点数\n",
    "        color_fit = double_sine_model(phase_fit, *popt)\n",
    "\n",
    "        # 优化图像清晰度和细节显示\n",
    "        plt.figure(figsize=(10, 6), dpi=150)\n",
    "        plt.scatter(phase, color, alpha=0.7, label='Observed (clipped)', color='darkcyan', s=30)\n",
    "        plt.plot(phase_fit, color_fit, 'k-', lw=2, label='Double sine fit')\n",
    "        plt.xlabel(\"Phase (folded by 0.394 days)\", fontsize=12)\n",
    "        plt.ylabel(color_label, fontsize=12)\n",
    "        plt.title(f\"{color_label} Color Variation with Phase (3σ Clipped)\", fontsize=14)\n",
    "        plt.grid(True)\n",
    "        plt.legend()\n",
    "        plt.tight_layout()\n",
    "#         plt.show()\n",
    "        plt.savefig(f\"{color_label}.png\")\n",
    "        plt.close()\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"[{color_label}] 拟合失败: {e}\")\n",
    "\n",
    "# 重新绘制所有颜色组合（含 3σ 剔除，增加图像清晰度）\n",
    "for b1, b2, label in color_combinations:\n",
    "    analyze_color_sigma_clipped_v2(df_all, b1, b2, label)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0f80229d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import plotly.express as px\n",
    "from datetime import datetime\n",
    "\n",
    "# === 设置真实周期（单位：天） ===\n",
    "true_period = 0.394\n",
    "\n",
    "# === 加载 CSV 数据文件 ===\n",
    "file_path = \"/mnt/7b21f1e1-eb25-4cd5-bdb5-06d7d82fa253/Temp/force_photmetry/result50cm/asteroid_photometry_results.csv\"\n",
    "df = pd.read_csv(file_path)\n",
    "\n",
    "# === 函数：准备指定波段的数据 ===\n",
    "def prepare_band(df, band):\n",
    "    df_band = df[(df['band'] == band) & (df['calmag'].notna())].copy()\n",
    "    df_band['timestamp'] = pd.to_datetime(df_band['obs_time'])\n",
    "    df_band['mjd'] = (df_band['timestamp'] - df_band['timestamp'].min()).dt.total_seconds() / 86400.0\n",
    "    df_band['phase'] = (df_band['mjd'] / true_period) % 1\n",
    "    df_band['tooltip'] = df_band.apply(lambda row: f\"time: {row['obs_time']}<br>mag: {row['calmag']:.3f}±{row['magerr']:.3f}\", axis=1)\n",
    "    return df_band\n",
    "\n",
    "# === 处理 g 和 r 波段 ===\n",
    "df_r = prepare_band(df, 'r')\n",
    "df_g = prepare_band(df, 'g')\n",
    "\n",
    "# === r 波段相位图 ===\n",
    "fig_r = px.scatter(\n",
    "    df_r, x='phase', y='calmag', error_y='magerr',\n",
    "    hover_name='tooltip', title=f\"r-band Folded Light Curve (P = {true_period} days)\",\n",
    "    labels={'calmag': 'Magnitude', 'phase': 'Phase'},\n",
    "    color_discrete_sequence=['red']\n",
    ")\n",
    "fig_r.update_yaxes(autorange=\"reversed\")\n",
    "fig_r.update_traces(marker=dict(size=6, opacity=0.8))\n",
    "fig_r.write_html(\"lightcurve_r.html\")\n",
    "\n",
    "# === g 波段相位图 ===\n",
    "fig_g = px.scatter(\n",
    "    df_g, x='phase', y='calmag', error_y='magerr',\n",
    "    hover_name='tooltip', title=f\"g-band Folded Light Curve (P = {true_period} days)\",\n",
    "    labels={'calmag': 'Magnitude', 'phase': 'Phase'},\n",
    "    color_discrete_sequence=['green']\n",
    ")\n",
    "fig_g.update_yaxes(autorange=\"reversed\")\n",
    "fig_g.update_traces(marker=dict(size=6, opacity=0.8))\n",
    "fig_g.write_html(\"lightcurve_g.html\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "fcfdb41b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEYCAYAAACju6QJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABiMUlEQVR4nO29eZhcV3Xo+1s19VA9d6u7rbHlCdsYGbBkHMdE7Zj4MsRMQQ8njolJwCFECe8GSLjce8E4L1wSLlycmMAVhAAOiRwFBzBxgsBJi8EoyAZbYBvbst2SLLm71ep5rGm9P05Vqbq65rm61+/7+pPq1D7n7LNrn732WnvttURVMQzDMAyjvnBVuwKGYRiGYeSPCXDDMAzDqENMgBuGYRhGHWIC3DAMwzDqEBPghmEYhlGHmAA3DMMwjDrEBLhhFIGIDIrI89WuRzpE5DERGSzw3FeKyJOlrVHp7i8iAyKiIuIpwb1uF5G/K/Y6hlFJTIAbBiAivyEiD4nInIi8ICL/KiLXVrgO/1dEvpzi+A4RWRaRrnyvqaovVtWhHO+vInJhwrnfU9UX5XvPUpF8fxEZFpFXFXo9EfmiiASiv/GEiHxbRC4pTW0No/KYADfWPSLyR8CngI8CfcBW4K+BN5T5vu6kQ18E3iwi/qTjbwO+qaoTeVy7aK10jfIXqtoCbAbGcNrcMOoSE+DGukZE2oE7gN9X1XtVdV5Vg6p6n6q+P1qmQUQ+JSKno3+fEpGGNNe7VESGRGQqar5+fcJ3XxSRz4jI/SIyD1yXeK6q/hA4Bfxawjlu4DeAL4nIBSLy7yJyVkTGReQrItKRUHZYRP5ERI4C8yLiSdRaReQqEflhtG4viMhdIuKLfvfd6GUejWqob01eHsjh2T4tIv8iIrMi8p8ickGaNvqSiLw3+v9NUc3/3dHPF0a1Y0m8v4jcjTOxui9avz9OuOTNInIi2ib/PdU9k1HVBeDvgcsTDvtE5MvR+j8mIjsT6vwBEXkm+t3jIvKmhO8uFJFDIjIdrcM9Cd9dEtX0J0TkSRH5f3Kpn2HkgglwY73zC0Aj8M8Zyvx34GrgpcAVwFXA/0guJCJe4D7gINAL/AHwFRFJNEP/BvBnQCvw/RT3+jKOxh3jVYAX+FdAgP8FbAQuBbYAtyed/+vA64AOVQ0lfRcG/ivQg/Pc1wPvBlDVX4qWuUJVW1T1nsQTc3y2Xwc+AnQCx6LPmYpDwGD0/7uBZ6P/AvwS8D1NivGsqrcAJ4Abo/X7i4SvrwVeFH2eD4nIpWnum/g8LcDNwE8SDr8e2A90AN8A7kr47hnglUB79Bn/TkTOi373pzjt0omj2f9V9B5+4Ns4E4VenPb5axF5cbb6GUYumAA31jvdwHgKYZfIzcAdqjqmqmdwBvBbUpS7GmgBPqaqAVX9d+CbOAN3jK+r6g9UNaKqSymucTewW0Q2Rz+/Dfj7qFXgmKp+W1WXo/X4JOcEX4y/VNWTqrqYfGFVfVhVD6tqSFWHgf+b4vx05PJs96rqj6Jt+RWcCU8qDgGvFBEXjsD+C+AXo9/tjn6fDx9R1UVVfRR4FGeSlY73icgUzgSjBbg14bvvq+r9qhrG+R3i11HVA6p6Ovq73QM8jTORAwgC24CNqrqkqrGJ2a8Cw6r6t9E2/zHwVeAteT6fYaTEBLix3jkL9GRZM94IHE/4fDx6LFW5k6oaSSq7KeHzyUyVUdUTwHeB34xqiW8EvgQgIr0isl9ETonIDPB3ONp0ImmvLyIXi8g3RWQkev5HU5yfjlyebSTh/ws4AnIVqvoMMIcj4F+JMxE4HdXmCxHgOd03yv9W1Q5V7VfV10frku46jbF+ISJvE5FHossHUzim91jb/TGOdeRHUdP7b0ePbwNeETsnet7NQH+ez2cYKTEBbqx3fggs4QjKdJzGGYxjbI0eS1VuS1SzTCx7KuFzLun/voSjef8a8FxUcwPHfK7ADlVtA34TR3Akkun6nwF+DlwUPf+DKc5PRy7Plg+HcDRRn6qein5+G44Z+pE051QldaKIbAM+B+wFulW1A/gZ0bZT1RFVfaeqbgR+F8dMfiHOZOpQdMIQ+2tR1d+rxnMYaw8T4Ma6RlWngQ8BnxaRN4pIs4h4ReQ1IhJbZ/0H4H+IyAYR6YmWT7Vn+D+BeeCPo9cYBG7EWVfNh6/irG9/hKj2HaUVR3OdEpFNwPvzvG4rMAPMibN9KlmQjALnpzm3VM8W4xCOQIw5zw3hrKt/P2rCTkWm+pUTP87k4QyAiLydBOc3EdmTsOQxGS0bxrEsXCwit0TbzCsiu3JZozeMXDABbqx7VPWTwB/hOKadwdGc9gJfixb5/4CHgKPAT4EfR48lXyeA4wj1GmAcZyva21T153nWZ55zQvwrCV99BHg5MA38C3BvPtcF3ofjRDeLo1Hek/T97Tje7lPJ3tKlerYEDuFMKGIC/PtAc8LnVPwvnInUlIi8r8D75o2qPg58AsdaMwq8BPhBQpFdwH+KyByO89t7VPU5VZ0FbgBuwrFgjAB/DqTcwWAY+SJJzp6GYRiGYdQBpoEbhmEYRh1iAtwwDMMw6hAT4IZhGIZRh5gANwzDMIw6ZE0mPOjp6dGBgYGSXGt+fh6/Pzm3hJEr1n7FYe1XONZ2xWHtVzilbruHH354XFU3JB9fkwJ8YGCAhx56qCTXGhoaYnBwsCTXWo9Y+xWHtV/hWNsVh7Vf4ZS67UTkeKrjVTWhi8iroxl6jonIB1J83y4i94nIo9EQhW+vRj0NwzAMo9aomgAXJ03ip3ECQ1wG/LqIXJZU7PeBx1X1CpzsRZ+QaPpDwzAMw1jPVFMDvwo4pqrPRqM87QfekFRGgVYREZwEBRNApqxRhmEYhrEuqFokNhF5C/BqVX1H9PMtwCtUdW9CmVac0ISX4IRdfKuq/kua690G3AbQ19d35f79hYZoXsnc3BwtLZmSGxmZsPYrDmu/wrG2Kw5rv8Ipddtdd911D6vqzuTj1XRiS5UFKXk28V9wMhP9MnAB8G0R+Z6qzqw6UXUfsA9g586dWioHAnPkKA5rv+Kw9isca7visPYrnEq1XTUF+PM4yRpibGZ1isa3Ax9Tx0xwTESew9HGf1SZKhpGfTF8aJgjdx1h8rlJOrd3smvvLgZ2D1S7WoZhlIFqroEfAS4Ske1Rx7SbcMzliZwArgcQkT7gRcCzFa2lYdQJw4eGOfjeg8yPzdPS18L82DwH33uQ4UPD1a6aYRhloGoCXFVDOCkbvwU8Afyjqj4mIu8SkXdFi/0pcI2I/BR4APgTVR2vTo0No7Y5ctcRfH4fDW0NiEtoaGvA5/dx5K4j1a6aYRhloKqBXFT1fuD+pGOfTfj/aZx8uoZhZGHyuUla+lY6zvhafEw+N1mlGhmGUU7WZCQ2w1gvJK55z43MEQ6G6djaEf8+MBegc3tn9SpoGEbZMAFuGHVKbM3b5/fR0tdCOBhm4tgEAO2b2wnMBQjMBxjcO1jdihqGURZMgBtGDRLTrOWXhQOfPpDSmzxxzRuIa96LZxdxe910bu9kcO+geaEbdYvtqsiMpRM1jBoj0Zvc5XWl9SaffG4SX8vKyMLtm9tp6W/htoduY8+BPTbYGXWL7arIjmnghlFjJGvWsX+P3HVkhUDu3N7J/Nh8/HuwNW+j/kjUsn1+H4jTj+dG5mjqbsr6HqxnTAM3jBojlWadypt8195dBOYDLM8soxFleWaZwHyAXXt3VbK6hlEwK6xNHhenf3ya0w+fxuV1sTS9xPTwNAvjC/HytqtiJSbADaPG6NzeSWAusOJYKs16YPcAN3ziBvy9fuZG5/D3+rnhEzeYdmLUDYnWpukT03h8Hjw+D9PHp2lodTTuqeNT8fJmYVqJmdANo8bYtXcXB997EIAGGuKadSpv8oHdAyawjbolMXZBcCGI2+dGRAguBNlw2QbGHhsjMBtAI2q7KlJgGrhh1BiJmnUkGDHN2lizJFqbvM1eNKxEQhG8zV6au5vpGOigob3BLExpMA3cMGqQmGY9NDTEwGUDHLnrCAffezDlVhrbamPUK4nWpvat7Yw9NgZA18VdLM8s4/a6efNX3mz9OQ2mgRtGDROYC2TcSmNbbYx6ZoW1KRRh48s3svHKjQVZnoYPDXNgzwH27dzHgT0H1sU7YBq4YdQw82PzGbeU5brlzDBqlVL4cSRHJYxNZNe6yd00cMOoYcLL4YxbynLdcmYYa5nkTHzhYJip41Psf+P+Na2NmwZuGDWMu8FNYC6QNliLBXMx1iuJvh+Tz0zSfUk3DTSwcHaBM4+fQVwCEeLa+I5bdnDy+yfjviJbrt2y4nM635Fa9jExDdwwahh/rz9jsBYL5mKsR5J9P1weF2M/G2Ph7AJTw1O43C5EBK/fS0NbA+FgmKHbh+Llx58a54EPPsD4U+MZfUdq3cfEBLhh1DC+Fl/GYC0WzMVYjySbzLsu7ALg7NNnCcw729Ii4Qgd2zoARwvXkMbLL5xZwO1xs3BmAXEJDW0N+Pw+jtx1JON90pWrFmZCN4waJ5uTjwVzMdYbiQFgAJp7mul9cS/jT47jcrkQl7Dh4g009zQDzrJSQ8u5ZabgQhB3g5vgQjB+LJXvSPJ90pWrFqaBG4ZhGHVFqnDDngYPF7zqAt76tbfSvq0dt88dX1Zye9w09zbHy3qbvYSXw3ibvfFjqXxHcg1rXC1MgBuGYRh1RSbfj1TLSrtv343b646Xb97QTDgUpnlDc0bfkVr3MTETumEYhlFWSu3JHRPSidcc3Du4wjck+fr9L+2Pl++5uIeX/fbLVnihJ56f632qjQlwwzAMo2yUK8hKvr4fqcpf/Z6rS36fSmIC3DAMwygb2aIF1vI+61rH1sANo4Ksx3jNxvomU7TAWt9nXetUVYCLyKtF5EkROSYiH0hTZlBEHhGRx0TkUKXraBilwgYrYz2SyZO7VvdZ18tEu2oCXETcwKeB1wCXAb8uIpcllekA/hp4vaq+GNhT6XoaRqmo1cHKMMpJJk/uWozlX08T7Wpq4FcBx1T1WVUNAPuBNySV+Q3gXlU9AaCqYxWuo2GUjFocrAyj3GSKFliL+6zraaItqlqdG4u8BXi1qr4j+vkW4BWqujehzKcAL/BioBW4U1W/nOZ6twG3AfT19V25f//+ktRzbm6OlpaW7AWNlFj7nWPy2UkiwQjilvgxDSsur4vO81MPWNZ+hWNtVxzlar/AXID5sXnCy2Enc1gojNvjRtyChhWNKG1b2lZNdivF+BPjuLyrddtIMELPpT05XaPUbXfdddc9rKo7k49X0wtdUhxLnk14gCuB64Em4IciclhVn1p1ouo+YB/Azp07dXBwsCSVHBoaolTXWo+slfYrhafssJzbTuNr8RGYCxCYD2TcTlPq9ltPHr9rpe9Vi3K03/ChYQ7evvIdmH1hltaNrXHNu9p98sCnDzAzNrMiw9/yzDIA2qs5vTuV6nvVNKE/D2xJ+LwZOJ2izL+p6ryqjgPfBa6oUP0MAyjdmli1E4/U09qeUf+kcgRLZZ5uPa8V/wY/tz10G3sO7Kn6hDLVmv3sC7PMjszW3LtTTQ38CHCRiGwHTgE34ax5J/J14C4R8QA+4BXA/6loLY11T7Z9rPlQzaAQ+T7HetLWjdKSLnjL8swy3Rd1ryhba34gqaKvIYBSkjGglFRNgKtqSET2At8C3MAXVPUxEXlX9PvPquoTIvJvwFEgAnxeVX9WrTob65Naz0iUK/k8R7miZxnrg3STxbmROSczWIJ5utpOa6lInmjv27mvJseAqkZiU9X7gfuTjn026fPHgY9Xsl6GkUjn9k7mx+ZrftDJRj7PUUqrg7H+SDdZ9Pl98XzdiX4gg3sHq1DL7MSsUBPPTDB1fIrui7pp7j6XorTaY4BFYjOMLNR6RqJcyec5bMubUQzptof17eirqh9IPiT6jPS8qIfQUojRn44y8ewEJw+f5NRDp5gfm6/qOrgJcMPIQrWdz0pFPs9Ri/tzjcoQcz4bf2K84Chk2dJ97jmwp2ac1tKRaIXyb/DTd3kfLreLs0+eRRD6XtIHUFVnNktmYhg5UIzzWS05g+X6HLv27uLgew8C9WHqNEpDou9Dm7eNmbGZgnwfaj0NZy4kLwM09zTjafKgYWXz1ZtXlK3W0pIJcMMoI/XqDLYWBmAjfwrxfUg3Qa3lNJy5kM5npKGlYUW5ai4tmQA3jDJSz85g9T4AG/mT746Lep2g5kIqK5Tb46a5t3lFuWouLdkauGHkST6ZiswZzKgn8vV9qKe44fmSymdk9+27cXvdNePQahq4YeTI8KFhhj48xKkjp/A1++i8oDOrxrFWtqAZ64NErbOBhriASuf7UO4YCdX2H0llhep/aX/NLC2ZBm4YORAzFY49NoanwYNGlPGfjxMOhDNqHGtlC5qxPkjUOiPBSNYdF+XcrVCroX9ryYveBLhh5EDMVBgJRnB73bg8LlxuF1PHpzJqHKXagpaP2d4wiiEmoHou7ckqoMo5QV3L5vlSYSZ0w8iBmKnQ2+wlvBzG5XEhbiG4EMyqcRTrDFYNR6Fqmy6N+qCcuxXWSgjjcmIC3DByILaW3THQwZnHz8SPu73usu+PLtSTvVAhnGrCcN9t99Ha30pgvjZSPhq1Q7l2K9ST/0jyu9b+W+0Vua+Z0A0jB2KmQrfXTc+lPeCC4FKQDZdtKPuWmUI82YtZP0w2XYaDYWZPz3Lm8TM1tRZprG3qxX8k1bs2c3KmIu+HCXDDyIFk557tg9u5+f6bufXQrWXXRAtxFCpm/TB5wjA1PIXb6yYcDNtaZB2wVvwl6iWEcap3TVxSkffDTOiGkSPVCmxSSFjTYtYPk02XwYUgLrcLT/O54cLWImuTXP0lilleqaRvRD0EE0r1rolbKvJ+mAaeA2tlRmvUJ5k0kXR9s5jtPcmmS5fXRSgQomNbR97XMipLLpaXmJAff2qcmVMzPP1vT3PPG+/h8J2HM167Vrd1VRtfi4/nf/Q8x793nNMPn2bh7AIa1oq8HybAsxCYC1inXafU0sQt1d7TTANqMeuHyROG3hf30npeK26fu6bXIo3c/CWO3HWEcDDM9PFpIoEIvmYfGlGGbh/K2MdtW9dqhg8NM3t6ltBSCJfbFU85Gg6FK/J+mADPwvzYvHXadUg9aBuZBtRi1w8TJwy3Dt3KjZ+7sebXIo3cLC+Tz00yPzaPy+3C5XFEgLvBjYY047hWTFjgWpoMl5Ijdx2h9bxW+i7vc9owongbvbg97oq8H7YGnoXwcthiWa9D6iEJSbZ17lKuH9bDWqSRm79E5/ZOxp8cx9d8blzTsGYd1wrd1lXqOAaFrMOXa+0+9g6KS2jucZKcaETRiBZ97VwwDTwL7gZ32UIFGrVLPSQhKWcYS6P+iAmp5ZllJp6ZYOLYBOD024PvPRjXfHft3YXb4ya0FEJViYQiRMIRmnubM/adQpdlSml6L8QyVk5rWrp30N3gLvrauWACPAv+Xn9d7EU0Sks1hWOiuXHy2cm0A0297JM1yk+ikOq+qJuuC7pQVWZHZkFZIbgAdt++G1wQWgzh9rlp39aO2+vO2HcKXZYp5WS4kMlAOdfu072D/l5/0dfOBTOhZ8HX4itbqECjdilk61YpSDY3RoKRtObGbGEsqxUO1cKwVp5USz5nHj+DovRc3BM/Fiu758CeVVm1cvmdMi2lpPvdSxlRrZDtkeUMyZruHRzW4aKvnQsmwHPA1v/WH+WM8ZyJ5IFY3LLCMS1VPZOPDx8aZuhDQ5x66BTeZi9dF3ZVJH567N6VjttupBZS4UAYZeVabCYfiZjlp5BwoJl+91ST4dmRWQD27dyX1ySvkMlAuUOypnwHh4ZLcu1smAndMNJQjbSBxZobYwPpmcfP4G30QgTGnxgnHMyc9rRU2Faj6pBqycftc+P2rVyLTSe4ig0Hms+OCARi84p816QLWTZay0tNVRXgIvJqEXlSRI6JyAcylNslImEReUsl62cYlabYtffYQBoOhhG3nEt7Opw57WmpqAfnv7VIKiHV2NlIU1dTToKr2HCg2X73xMmwf4Of1vNaC5rkFbIOv2oCwWrHvnqlaiZ0EXEDnwZ+BXgeOCIi31DVx1OU+3PgW5WvpbGWqNW12cR6+Vp8zL7gmBd9LT40rHmtvSenPRWPI8RzSXtaCuopg9RaIuWSzx2DADktAyWa4BfGF5g6PsWm6zcx/J1hhg8NZ31P8vndi12TLmRJM3bOWlviqeYa+FXAMVV9FkBE9gNvAB5PKvcHwFeB+rd3GBUnJhxHjo4wPzZP+9Z22je318yLmzygTD8/zewLs8ycnsHj9XDhGy7Mq47xtKfbOjjzhJP2VNUJh1oJJ7xqOf8Z6QVbPmvL4UCYM0+cweV2jLMujyun9ySf372ak7x6iO+QD6JamQ3nq27smMNfrarviH6+BXiFqu5NKLMJ+Hvgl4G/Ab6pqv+U5nq3AbcB9PX1Xbl///6S1HNubo6WlpbsBY2UVLP9AnMBZk7OIC4htByKB1fwNHpweVyEA2EioQhurxt3gxt/r3+VGbDcTD47SSQYQdxCJBQhtBQCQFyCp8GDq8uFv+lcvQJzAWegXQ6nrHPiM6sq4eWwEx3K76V1Y2tFni9bHSuFvbu5E+s34UA4fsx3no/wRBgRweV10Xl+ZgGb6++e2EfFLWjYCXzStqWt7P1k/IlxXN7VK8eRYMRJE1wiSt33rrvuuodVdWfy8Wpq4JLiWPJs4lPAn6hqWCRV8YQTVfcB+wB27typg4ODJagiDA0NUaprrUeq2X4H9hyIz/SPf+84ngaPI7Ab3HQMdDD22BhEYNsvbSMwF+DM/JmKa+T73rcvHsnp9MOnHbO3WwgHwmx75TZ8v+pj+jvT7Dmwx9HWb3e09ZiWk6rOtbpUUGns3c2P4UPD3PPGe4hEIvj8PrZ9aBvBfwmiEWVudI43PfSmlOfUQ1azGAc+fYCZsZkV2v/yzDL+Xj+DvzdYsvtUqu9VU4A/D2xJ+LwZOJ1UZiewPyq8e4DXikhIVb9WkRqmwAbH+iFxrS2+JuwWggtBpoanEARvqzfuSAOVN6UlmhODC8H4JMPb7AVWpiXM1fyX6xqh9WUjkYHdA5z/qvPj/TEWJz2b53oh68n1lJq3lqmmF/oR4CIR2S4iPuAm4BuJBVR1u6oOqOoA8E/Au6stvGs9wYVxjsQ0f+HlMKHlEOHlMN4mL8uzywArUmRWw1s60XvY0+QhtBQiEo7E65WYlrCUHt7Wl41UJPZHIG/P9VrfMlhskp9ao2oauKqGRGQvjne5G/iCqj4mIu+Kfv/ZatUtHWvNAWKtMnxomKEPD3HyhyfRsOJt9qIuZ50tIhF8rT7ELTR1N8UTEEB1vKUTvYfnRuaYW5qjfWt7fPuPN+KND56ldP6xvmykIrE/RoIRx7Scg+d6jHrYMriWAnNVNRKbqt4P3J90LKXgVtVbK1GnTNRrh11PxDTLqeNT+Pw+IsGIE++5wY3P72PDZRu49dCt8XLLM8tVNaUlmrH7dvSx89qdnPz+ybhZu21LW3ywKaX5z/qykY6YgBsaGsq4Llwpb3Jb6kmPRWLLA8v+VPvENMtI0PEu9zZ7aWhroKGtgU1XbSIw7/x+tWBKS2XGPnr3UXbt3RWP/pZoMi9lna0vG8Wy5dotjB4d5bmh5zj10CnOPHGG0aOjjDw6UrIAKbbUkxmLhZ4Ha80BYi2SHMjE5XHFHdeSBVS1TWmFmLFLVWfry/XJ4TsP8+DHH2RxcpGmziauef81XP2eq9OWL5f2OnxomKN3H6VtaxsLYwssTS+xMLZA20Ab3Rd1lyzOgi31ZMY08DyoBa3NyExMs+wY6CASjhAJRdCw4va6ay7+caXDjiamKT1y1xF23LLD+nIdcfjOwzzwwQcIzAbw+X0EZgM88MEHOHzn4ZTly6m9xgRrx9YONu7cSGNHI74WH+HFcEkd2iw0b2ZMA8+TamttRmZimqXP76Pn0h4mjk0QXAiyaecmBu+orTSwlYxIlWrLz9G7j5rQriMe/PiDuD1uPI3OsO1p9MCSczyVFl5O7XXyuUlcHhdnHz5LcCFIcD6I1+8luBCMlymFoLXQvJkxAW7UHMWY/ZJjQm8f3F6zTi+VNGOXczDP5/cyh6TCWZxcxOdfqY26G9wsTi6mLF9OR0Wf38fpH5/G4/Pg9rkJLgRZnl2mqaMpXqYUgtaWejJjAtyoKQoJDpFKKOw5sKfoepRb0FQy53i5BvPhQ8Pcd9t9LE4sEg6EmTo+xcjREW7cd2PKPOVrKZFEIRTSr2LnRIIRFicXaWhtAIXgYpBwMIzH50mZcKSs2mtCYEwRJ+xvYCHgLFlFtGSCtpLvSD1iAtyoKfLVFMshFIYPDXPfO+9jaXKJcCDMxLEJnv7Xp2nd2Er/Ff0lFeaVWpIp12A+9OEhZk/P4vF58DZ6iYQizJ6eZejDQ9w6dOuKsql+28WpRe69+V5a+lvWvEZe6OQ0dk7nBZ2cfeosS1NLgBMvXxD8G/3cd9t9tPa3EpgPxNuxnNprYC5A7+W9TB+fJrgQxNfqo3VLK3MvzDE3OldSQVvMO7LWLT4mwI2KkcvLlK+mWA7T8NCHhph9wRFKuCAwGwB10izWq9ZYrsF89NFR3F53POymy+PCrW5GHx1dVTb5t10YX2B6eJpIJEL/Ff3xtt1xy44Ve+HXyqBbSF9NPKehrQERYfzn46BOW3de0EljeyOjPx1laWKJTVdtWtFHE7XXWAz9g+89WHTc8tiEcOOVG+PllmeW6XtxX9HWr1KRz4SpXgW9eaEbFSFXj9h89yeXyks10UP75A9PIjh5tEOLIcQlcUFeD+EiU1GuHRSKkpxoSETQVXmJVv+2U8enAGhoPReKMxwMc+j2Q2ty328hfTX5nK4LumhobcDX5uOCX7mArvO7mBqewu11Ew6u9gAf2D3AngN7uOETNzgxEJS82jXde7vl2i3xkKsa0YwhVwsh8X0sdE95rqFe63mvuQlwoyLk+jIlxmLOZWAoRUCS5BdYI0pwMeikIA2vFlCpBt1SDDjlJjaYx4LElELD6N/RTygQctY+VZ2UqIEQ/Tv6V5VN/m0DswEUpWOgI15mYWyBcChcV/G1c6WQvprqHLfPjdvnjn8OLgQRkXgCHFjdRwuNW57uvJPfP1m2LbWlEqi5TpjqMaZ7DBPgRkXI9WXKV1PMV+CnIvkFbmxvRFWdgdEtTh5xPWfyTB5063kGXyyDdwzSel4r4hIn25tLaD2vlcE7BleVTf5tG9ob6BjooLl7ZTz6tbrvt5C+muqcxs7GeKx8jSgur4tQILQiMU9yHy3UUpXpvHJMCKF0AjXXCVM97zW3NfAqEku6MfroKIrSv6O/5vYql4p8nKjycVophZdqbG124ewCU8NThJZDgJMJzOV1OTnEvW66LuqKD7qXXHsJB/YcYPK5SeZG5mjqblqX0aIGdg9w4+duzHn9MPG3TRWPXjyCv9e/4py1su+3kL6a8pzo5Ch2rPfFvcyensXtc6f1AC/UiTGW0S8SdFLcdgx04Pa6y/p7lGrHRK5+H/W81zyjABcRF/AWVf3HCtVn3RDbfjN7etZxAhIXp398mvveeR83fm71Fpx6p5wescV6cndu72T8qXGmj0/jcrvw+X2OOTgYoX1zOy3ntYBCYD6Af4ufS669hKN3H407x4w/Oe5Ex2r2xbOb1csMvhQU2v4phdPtgxy9+2jVk8yUi0LaKnZOzNEqlRNashNW8sSgkPcvMBdg9vQsoaUQHp+T6nb0p6O0bkxtYSkVpRKouU6Y6nmveUYBrqqRaMpPE+Al5shdR1icWMTj88Q9eD14WJpcWpOaWy3v59y1dxf3vPEeUBC3nNO4L+yi5+KeVV61B/YcWOlN3NpAcD7I1PGpuACvlxl8tUkl0Ppf2l+T/aSaZPOozjYxGNg9wI5bdqyKo57pnPmxeVrPa6W5u5mp41MEF4J4G7209reW9fcopUDNZcJUy2NTNnIxoX9bRN4H3APMxw6q6kTZalVDBOYCcVNpKbcXTD43STgQxtt4zvFE3M464lrV3Go1DO3A7gGaNzQTmA0QXAzibfbSva2bpq6m+G+RqOFMPDNBz4t6aMAR4B0DHYw9NuY4ZZUwiMV6pVb7STUpdrtkLPlI1wVdcaF49O6j9L+0P+354eUwvhYf4pL4xFQjytzoXEmeKR3VEKj12udyEeC/Hf339xOOKXB+6atTWwwfGmbm5Mwq56R8vS3T7aOcOj5FJBSJa+AaVty+8q4vGamJ7UNONNstzyzTub1zlfYzdXyKscfG6Lu8j+aeZpq7m+kY6GDx7GLJg1gYhVOuyXc5Sbcfudh14UKC6Lgb3ATmAlVZG65XgVppsnqhq+r2FH9rXniD0+ljXpCFekNm2kfZ1NVEKBAiHAzHt980djbWVMasYhg+NMzks5M1u7UqcevX/Ng8syOzKT2Ek71iuy/qBmDi2ES8rNvr5s1feXPJPXKNwkg3+a61PphIpt0MxW6XTPa0jjlsLk8vp20ff6+/rHu9jeLJKsBFxCsifygi/xT92ysi3mznrQUmn5tE3On3AGfb+zt8aJh7b76Xs0+e5ezTZ+PxoqePT3Po9kO09rfSdUEXKEQiETa+fOOacWCLDUaRYKQmB9DkwRJw7ErCqu1ryYNfc3czvZf3EglFLBVnjVKKyXelybR9qtjtkquC6AxPIQi+Vl/a9vG1+KqSPrkeYirUCrmY0D8DeIG/jn6+JXrsHeWqVK3Qub0TDa+MKBWb9WZzKol9vzS9hLfJS3g5zOhRZ7uYt9FLJBIBnHCIb/3aW9fc4B8bjMQt8QEidjxdTPNE0+GWa7dw8vsnGTk64qQqbPaWNA75kbuOEA6GOfu0kw7R5XYRDoWZPzPPBa+6YMV9UnnFenwezn/V+TUTNnItU0iYy8nnJulx96w4Vis7Awoxkxe7LpzsGLY8u4zL5VqxdzxdXIZKjk2W8CY/cgnksktVf0tV/z3693ZgXdhQdu3dFZ/tZjOrJs9g49+3Njj7iT0ux1QedAR37Nxa1woKJZ/gCMna8PhT4zzwwQd44ScvMDcyR2A2wPzoPONPjZdMix95dITp4WnCy2EQWJpeIjgfJLwcXmUtKEWwGKMwCg2SE5t8L4wvcPrh0xz/3nFO/ejUqnSclaYYM3kxgVOSg+g0tjfSPtAed05Lvlchz1WM1hw7f/8b9zN1fCplWFhjNbkI8LCIXBD7ICLnA+HyVal2GNg9QNuWtpQmpGwCKvZ9x0AHkXAknmZPI0okHImHj4ydU6oXoFpmp+T7+/y+nNfskidDC2cWcHvczJyawe1x42l0ttotnFko2cscXAgCrIx3Lo6XbfKgUa444kZ2Co3KtWvvLsLBMKM/G41HiAsuBZkdma2qSbacZvJsJE4A3vyVN+P2uktyr2IjESaeT8Rx5j3z+BkWzi4AtWM5SabaYy7kZkJ/P/AfIvIsThbYbcDby1qrGsLX4ktpJk1lVp1+fprFs4vs27mPuZE5wsEwHVs72HDZBqaGpwBnq9iGyzbEw0fGQkcWYzaqptnp8J2H+e6ffpfFiUVcXhddF3bFHcJQaA23Zt1alWw6DC4E4x6wrjZnjilucdIWluhl9vq9LM8uEwlFiIQjTrzzaIYncH73kUdH6s6Lea1RqPf1wO4BHn7yYTyNnngUsZ5LenB73Tkv45Tj9y6XmTzfuifea+TREWePt9+7YtKaK8VucUs83+t3lhvFLUwNT9Hc3VyTMRVqxdSfixf6A8BFwB9G/16kqv9R7orVOsmz5akTU0wcm6Cpu4mWvhaaupuYODbB1Ikpmjqb6L6om84LOunY5oQiTJz1ohSkZcSoVjD+w3ce5oEPPsDy9LLj7BeBsz8/y8LZBVr7W2nd2IrL68qqtSabDr3Nzkvs8jlhTMGZlXubvSV7mft39DthIRvcCBK/b2zwmX5+moUzC3XlxbwWKcb7WsPK5qs2s+2V29h45Uaau5tzXsYp1+9dDjN5oXUf2D3Arr27aGhroOuCLrov7C7ouYuNJZ54fsc2x2IJxMfXWlyuqpUEKGkFuIi8OfYHvA64ELgAeF302Lom2ay6eHaRrgu76NjagbiEjq0ddF3YFd8b7O/1c+O+G7nxczeCwMkHHQctn9/H7MhsyV6AQs4vlAc//iBuj5MVScRxVhOXMPHMRDxYROf5nVkHo+TJUPOGZsKhMG2b2giHwoSWnGxXzRuai36ZY2avkaMjTJ+YpnlDM+ftPA93gxtVpX1rO8szy0yfmKZta1vVX9D1zPChYebH5jn10CmOf/c4J35wgueGnmP06Chbrt2S9fyYFSeRXJdxsv3ehZpPy2EmL0aYlEIQFbvFLfH85p5mNly6AXEJLperZperaiUBSiYN/MYMf79aipuLyKtF5EkROSYiH0jx/c0icjT696CIXFGK+5aKxNlyS38L7ZvbV3zfvrmdlv6WVQIsMBegb0cfW6/ZCjghC6efn15xbqEvQCHnF8ri5KKjwbod8zMALogEInndP3ky1HNxD9d/9HrOe9l5tPa34mv14e/z03NxT8aXOZdtfTFNpfvCbtq3tjNzYobFs4tsvHIjG1++kUgogr/Xj7/Xv+r3rNW1uLVI7LcCaN/aTmA+wNL0Em6fm7atbRy9+2hWoZnPPuZinC7z0VrL4U9RjDAphSAqdlKSfL7b56Z9Wztv/dpbazamQrXG3GTSroFHvc3Lhoi4gU8DvwI8DxwRkW+o6uMJxZ4DdqvqpIi8BtgHvKKc9SqUXAPwp1ovigmSpo6mgmL/5hI7uBzre02dTQRmA3gaPc5SQARUFfFI/P7DOpzTtVJtV7n6PVevqnu6NbrYoLowucDcqTle+MkL/PzrP+flv/tyXvdXrwNWt33H1g6aOprw9/pTxjuv1wxFa4HE3+rs02dpbG8EHK26Y2sHyzPLWddYY/uY060pJ/arRJ8VgIXxBSaOTRAJRTiw58CK96XYNd9Sb80qJvlHPuemG0OK3eJWj7HIayUBSlYnNhFpBz4M/FL00CHgDlWdTn9WTlwFHFPVZ6P32Q+8AYgLcFV9MKH8YWBzkfcsG7n+oKmcWNo3txNaDOHv9ZflBUjlcHHfbffR2t9KYD5QsEC/5v3X8MAHH8DtceP1ewnNh9CI0nt577n98EPDeV0zmVjWtlgQnKnjU4wcHeHGfSsD3hy56wgLkwvMDM+AC1xuFxpRHv7Mw3Rf2M3V77k6L4eoXH/PSjg+rUcSf6vgQhC3z42IxHcP5KolphOWsXciHAyzMLbA0vQScy/MsTy3TMuGFsYeGwOg9/LeVQ5K2fpRpftEMcIk07mJk+9ik6lko95Cp9bKpENUNXMBka8CPwO+FD10C3CFqha1Di4ibwFerarviH6+BXiFqu5NU/59wCWx8im+vw24DaCvr+/K/fv3F1O9OHNzc7S0tGQviDNznR+bJ7wcxt3gxt/rX2Wemnx2kkgwsiLCWyzvdOf5hWt3me6dfM9IKBLfOuX1e9Gws72tbUtb/CXO9hwx5sfmmRuZQ8OKuIWW/pYV+Zzn5ubw4cv5esmcfeoswfmgs80rikYUr99L98Xd8WPjT4wTXAyuvoA66/OeJk98b2lsO5+4BJfHhafRk7Lts7VDYC7AzMkZZ+3fLavasRTk0//WEol9NrgQPLdEI46zYeyd8ff60/5Gmdpu8tlJQkshwoGwswMBp1+patyfw93gXpGnIPaOpnqHwwEnHLK4otnsfG4nP3eaPpHPO5aOxGvElrE0onlfL11dEtuvXOPWWqXU7+111133sKruTD6eyzayC1T11xI+f0REHilBnSTFsZSzCRG5Dvgd4Np0F1PVfTgmdnbu3KmDg4MlqCIMDQ2R67WGDw1z5EvnZmQX7714talXzs1kE2e8xayDDR8a5uDtK695Zv5M/Jr73rePlr6WuBA8/fBpQkuOtrztldsAJ3GH9qozI89wrXw5+M2DDN8+XPD1/vxNfw4Cbq87fiwcDIPCr02e65YHPn2Ap+59CpfbFR+QI+GIE0nPBRe9+iLGfj7G9HPTeJo8+Pw+wsthwqEw13/0eq4evDrvZ0tlZo+1Y6kitOXT/9YSie9JaDm0QiP2+Jwlmx237ODoXx9N27cytd2+9+1j9tQs4UA4LqQBAgsBXB4XW6/ZSsQViR+PZeF600NvWvUOTz8/zcSxCbou7GJhbMGZcAjx7aLJfSLb+5pT+6S4RrHjSDKJ7Zc8hiS3ibGSSr23uQRyWRSRuOAUkV8EFktw7+eBRFfSzcDp5EIisgP4PPAGVT1bgvuWhVwdW8rhxJLNkzTZ4SK4EERE8DafC2kfMwHm6pWaqxfu/Nh8UV6uisYFcgwRQRPmesOHhpk/4wSBiAQjhMNhVB3NJyb8xSWEF8N4mjxoRJ1Urn4vXRd2cfL7J3OqSzK14om6Fkl8TyIhJ0/Axis3EglG4u/Mye+fLLhvdW7vZGl6ieW5ZRYnFp1IfItOnAFBsm71SrcDJbjoxDBwuV3x2A+hQIhnv/Ns/F0Z+vBQ0Z7fld7GVCtOW8ZKctHA3wV8OboWDjAJ3FqCex8BLhKR7cAp4CbgNxILiMhW4F7gFlV9qgT3LBv5OLaUer0n25pc8jqXy+sitBSKm6ATHXaSc10nXwvyC2IQyymcrm7Z6N/Rz+kfn8aDJ26mDgVCbHz5xlV1advWxszxGQhDBEd7EhE6L3AGmeBC0NG8A+G45UEjWrDALcZ5yMhOtvfk4HsPFpxic8u1W3ji3idAiPer4HwQf6+fthe1xYVVujXlxLrt27kvXo94DAOPi+BCkIWzC4z9bAxvozf+rpx66BR9O/oyvmPZKDa9aL6Uw2nL/EeKJ5dALo+q6hXADmCHqr5MVR8t9saqGgL2At8CngD+UVUfE5F3ici7osU+BHQDfy0ij4jIQ8Xet1xUUxtLNTuefn6aqeEpPtb5Me554z1EQhEWpxY5+eBJQkshUMdcOH9mntGfjRJcCtJ9STcuj4uxx8ZYGF+IXytZKOWjpYeDYYa/O8zph0/HQyPmI+QG7xik9bxWR4OOhsRsPa+VwTsGV9Wl/yX9bHjxBlw+F6KCp8FD20AbXed3AecG10TLQzEC12KkZ6ec4SaL0QpPfv8kbQNtjjUnuszi8rmYPjHN7OlzcRnSWckSn2tuZC6+DTQWiCS0FMLT5OHs047RsOvCrnPvSrOPiWMTBdW7FM9eCKW2HFYqcE65qXY41bQauIj8ETCtqn8DoKoz0eN/ALhV9VPF3lxV7wfuTzr22YT/v4M6yXpWTW0seXY8/fw04z8fx+VyOQ4/qpx96iwaUbpf1E375namn59m+sQ0GlI8jR66L+qmubsZEWH0p6NOVLmuppVeqdEZ81P3P0VjWyMdAx3xZAjptPTO3+7E5XIRnA8y9tiYE/3M68555j6we4AbP3dj2pl6sibSud2Jdjc3OscNn7iBg+89yPLMMr4WH80bmpmYnqB9Q3vW8K651q0WPFEzUU0tp9zhJovRCiefm6T3kl5ae1uZOj7F8syyM0H0OPneM60pJz9XOBiOC+T2ze20b2tn+sQ0jW2NzI3N0fvi3hVJQzov6GT0p6PxfpnqHcv2e1VjG1MpLYfFbsWrBTL170qRSQP/beDuFMf3Rb8zEqimNpZqTc7b5MXb7MXlceH2uuMxvxfOLMQjxfXv6MflcbH5qs3x2Ozpcl0D8RlzQ1sDwYUgZ544E9fU02npbp+bDZdtwOv3QgQWzy7mPYBnCi+ZSRNJbpemzia6Luxi+sQ0Jx88CULRwqSYDFHlJh8tpxyaRLnXaYvRCmP9prmnmY1XbnTicDd7aepsylrX5OdKjrrYc3EPN33tJvY+uZcLXnUBwcVgPCPa6YdPE1wMsmnnJvy9fiaOTTDxzATLM8sMfXiI+955X06/V70n2FkL/iO1EE410xq4qmogxcFlSfYqMiqijWWanSevyU08PbFiywfqBFlZml5yBpGFIN4mL6FgiMBcgHAwzNTwlJMb2+ui74o+bh26NX56LMNYQ1sDnds7OfP4GVCYGp7C7XMTmA9wybWXxJN/xNbSwQmP2NzTHPdaTaXRFKolZtNEYu0SE2ZNHU20b253ys2t6t5rily1nHJpypVYpy1UK1zVb2YDhMNhXMsujn/vON5mL+1b21PWNV0sB7fXzW0P3bbi+JZrt/DU/U/h9rhxN7gJzgeZmJ7g+o9eT/9L+511/P4WfC0+nv/R84SWQo4lLCoQIL1WWm97pxNZC/4jmfr3BjZUpA4Z18BFpC+XY4ZDObWxfLSpzu2d8T2ocQRQJ8xpeDmM2+eOe6OfffYsoz8dJbQUwuV2HNxmT69Mu5g4Y27ubmbDZRvwNHtYmlnC3+t3tvTcfTRev9haeiwZCaR+QYtdC8tVE6mF2XKlyVXLKVfb1KrncmzCuDy7zMQzE5x9+qwTElidCa/b58RQP/3QacZ/Pr7KIpHPc538/km6Luxysmwl7XxIbvdIMILH52Hq+FT8/HrTSnNlLfiP1EL/ziTAPw78i4jsFpHW6N8gcB/wvytROeMcR+46QjgY5uzTZznxgxOcffos4WA45SC7a+8uGjsbCQWcJCCBhYCT4ScCEY3mJQ87W6w6z+8kshTB2+iNB4Hou7yP1vNaV1w7ubM2dzfTc3EPF7/2YvYc2LNqS0/3RY6He3g5nPEFLYXwyGXitBZMdvmS6wBTrrapxUE6OR5+1wVdNLQ10La5DfFE4wcEIwQXgqg670PypDKf55p8bpL2ze1svHJjPCta+2ZHs09u95i/SizaHNTGhKcc1PsSANRG/04rwFX1y8D/BO4AhnHikn8E+LCqfindeUZ5GDk6wtTwFOHlMJ4GD+Flx+Q9enR0VdmY49fGl28kFAgRWgzR0NqAu9EJR7k8sww4gSbaN7cTWAiw6apN59Iu9qxOu5itsyYPRrG1dFXN+IIWKzwO33mYT27+JH/m/zM+ufmTHL7zcMpytTBbrjS5DjDlaptaHKTTTRinj0/Te3lvPINZ7DsN66pJZT7Plaltk7/rGOggHAyvSjdcT1ppPtSy/0gu1EL/zrgPXFX/FfjXCtXFyEBwPogg8ahRLo+TKzswn3odd2D3ALceujUeLSwcCPPCT15wQotG18ZjUaKaOh1v80zrUdnW+FOtaQUXVoZATSRmxpx4ZoKp41NxL/hU905HLB+52+PG53fWMR/44APAuUQoMWol+UAlydUvo5xtU2vrtOnWLRXF4/Ow8cqNHP/ecTwNHickaoM7XiZxUpnLcyWmQ/U2O6bzWBS5WNsmtrvb66Z1Yyut/a3Mjc7V5K4GYyXV7t+5BHIxiqBU23i8zV4CswEn3nI08ETseKZ7P/udZwkuB9GgIt5ozGdVFicXmToxhdvr5pr3X8PRu48CmQfwTJ011Va2iWMT9Hh6Um6xiDlN9byoh7HHxhj96eiKMJm5CI9YPnJPo9ONPY0eWHKOJwvwdMIMiDvercVgErkMMPWwHa5UpHOe6t/RH58Me5o8hBZCINC9rTteJtukMvFd97X4mD09S+t5rfS9pI/JZyYZPTrKpl2bVmhpq9r9jrXZ7kZ5MAFeRkrp3dt/RT+nfnyKuVNzRIIRXF4XLZta6L+iP+O9XR4XOq9O6NGgMziFl8PgcrZ0vfkrb2Zg9wD9L+0vagBPFgKx8JJun5uQK7TCoxaImzEbaKDv8j4mjk1w9udnOf9V5+d878XJRXz+leZ3d4ObxcnUkX6ThVm59ynXE9XWJCpFOmtDbGJ55K4jzI3MMbc0R/vWdpq6muKm7EyTyuS+lOhR7t/gx7/Bz/LMMv4N/hXtXCvtnkrRMGqfjAJcRFzAW1T1HytUnzVFKYMVxLajeBo8uNvchJfDzI/Ms+XaLSnLx+7ddWEXp350yvF2EOKRyHou7SESdDzEEzXQYoRXuvCSMRLNkInfNfc009TVxNzoXF5JQBLzkccIL4dp6mzK6fy1EEzCyI9s1obYv8kCLdukMrkvJXqUpwt2VCukm8gO3D5Q7aoZWci2Bh4Rkb2ACfACKGXe4Nh2lIUzTrYjr99Lc3MzD378QY7efXTV+SOPjhCYDRBcDDomd9V4rrcNl21wsns1kJMGWsgyQMxU2Uhj/FiiGbIUe0Bj+chZcjTvWHaxa95/zYpy6epf6XjSa4l6jGOdXOdMk9XY8Vj5RAe2VCT3JW+zl9BSqKQe5eVq83QT2fmx+aKvvd6I/Ubyy8KBTx8o+3uRSzayb4vI+0Rki4h0xf7KVqM1RCYP1Hz3PydvR+kY6GD+zDzL08urzh8+NBwX9G6fG0+DByJOrOemnibcXifwCkrWLVyF7tOOeUDH8iEnetSWavvF1e+5mus/ej2+Vh+B+QC+Vp+TGjRh/TtT/dejZ3opqKc41rEIc3918V9xzxvvYfyp8ZzqnO8zltujvJA2zzW6XrqdIOHlcEF1XQsUEpkw8TdyeV0VeS9yEeC/Dfw+8F3g4ehfzSYVqSUyCap89z8nDxBTw1MIgq/Vt+r8I3cdoW1rGwhoWPE2e3E3uIkEIzS2Nca3OwTmA1m3cBW6TztmqnR5XSm3WPhafIweHeXEgyeA/EOaxl6wo3cfZcsvbOHm+2/mj57/o1XOa5nqXwv7OOuRegmKE5gLxAfUwGwAjSjTx6dZnFzMWud8nzG5L8Uc4kJLoYL7eDH1yTfwU6qJbMwDf71R6AQ1lSWj3O9FVic2Vd1etruvcTKtt+WbCjHZ+WZ5dhmXy0XHto6U57dvbsfn98XDoza0N9DY1sjeJ/fGy+cSzrAYM/PA7gGGdZg3PfSm+LHE9bYt12yJOxHlw+E7D3Po9kOEQ06q0nAwnNb5LFP915P3dSmpl6WHxFz0sTzdoaUQL/zkBTwNHjxNHuZG5lKem+8zJvalkUdHWDizQOf5nefC9ubYx0u13JOPf0c6xz5/rz+nOq81CvWNqcZ7YV7oZSadl2m+sYCThU1jeyNN3U0rshwlrzE3dzfH91YvzyyveiFz2f9b6pjFxTqODR8aZuj2IYiAr9lHJBBh+vg07dvaU14jW/1rxQu4nqjFONapBF9iLnpvs5fAXIDQYghVxdPmbBWbW5pj+NDwaue1ZyaZPj5N14Vd8Xcs2zPG+lIs9kK+fTzTroh82zwfYZJuIjusw2nrupYpVBBX473IxYRulJjEAA8nD59kfnw+J/NtYuSiN3/lzbi97pTm31xNw7lEEiq1mbnYyGtH7jqChjRu3nN5XLjcznpTqmuUsv7Vzv1bK9Ta0kM6k6e4JG4a7tjWERfe4hYnRr9A+9b2uIlzRZjVS7oJLgUZ/dmo42uSxzMW2sdLudyTr39HvUdFKyXJbbdwdoHnf/Q8E89MZHzvE38joCLvhQnwMpM86B++83Bc6+17SR+COOFQ80xtmUn45hPib2D3ALv27qJze2fc2zaxg6a61o5bdnDkriMZBVnsucefWJkMoljHsdjgmJioxeVxpU2UcuSuIyzPOEkrJo5NFBzusJ4ct8pNLYSQTCSd4EOID6hNXU1OrHMBt9fJDLbh0g3xuOTJ1/H3+Ol7SR+eRg/jT44XlKo0kVz6eCbBn2+b19okq55IbLv58fl4oqeeF/XknOI1EoxU5L0QVc1eKlZY5HZVvb1stSkRO3fu1IceKo2f3dDQEIODgwWdm2gSi5moR4+O0ra1jY6tHfFyMfN2PnugS0WqOsYCW6TqeLmUTyzTdksbM3fPrAiWkc/9kjmw5wDjT40zfXwal9uFuIXwchhxCW/92ltT1qGQ+6S6b7J5rBK/WzH9b70QizmQGLZXI4rvV33svGxn3DQ8NzJHU3dT2ncv3XXmRudWpQnNRHLfm35+mpkTMzRvaKb/iv60W4tK3ceK3Xa2nvterO2e/c6zuDyuFUspufwmpW47EXlYVXcmH89XA399ieqzLkilGYRDYRbGFlaUy9fRoZSm3Hy9W3Mpn8kbs1jtbdfeXbi9btq3tTtR3hZD4ILdt+9ecY1Se0qvx2xm9UImL+pcl50yXSffNczEPj5xbILpE9O0bW2j+6LujBpcqbVmM4sXTqztOi/oZNNVm1b4GtXSe5+vAE+dmcJISbpBv5hBotSm3HwFUy7ls5XJZrbPRGxw7Lm4h9ZNrVz46gu56Ws3rdo+VmqBa3vGa5d0gi/ZaTPb5LGUAjQmAPp29NG/o5+OrR1ZJ5LJ9UOcOA0H33twXftcVJNaf+/zFeBXlqUWa5RUP76/1494pOBBotSaZb4dNJfy2coUOwnJRbMo9Ytna4q1S7LGO/HMhLN+OTa/qk9l6jvlWNvPdyIZq98Nn7gh3n/Xu89FNan19z4vAa6qkXJVZC2S6sd3e90M3j5Y8CBRas0y3w6aS/ls3piVCARSDnNkLTluGSuJWXV8rT66Luii+6JuIsFI3kIvl8lhPktYhU4k6yVYzlqn1t972wdeRjIFCkk2+eZKqfca5hvMJJfyiWVi3piJZSoR8KAcQVpsz3htk+x7IW5Z4XtRCvLNYFdorvV6CZazHqjl994EeJnJ9cfP1WO00AGhFHXMpXzyc7RtaeOGAzesKFPOgAf1mGTDKA2VEHr5BiIqdCJZyDtifX/9kdWELiJ/lOLvd0TkpcXeXEReLSJPisgxEflAiu9FRP4y+v1REXl5sffMlXT7mMt1r1zXhGvZpJPqOWZOzqx6jnKtK9le7fVNtuRBpdi5UcgSVqJZPpYHIVs98n1HrO+vT3LRwHdG/+6Lfn4dcAR4l4gcUNW/KOTGIuIGPg38CvA8cEREvqGqjycUew1wUfTvFcBnov+WlRX7mL1tzIzNZDSTFUshs/paENiJDB8a5t6b72Xh7AIouNwuGtob8Kt/1XOUKwa55fde3yRbpzTsJBW55NpL8jJ7Z6IYzTgWI71taxvtm9sz1iPbO5Ksbc+fmbe+vw7JRYB3Ay9X1TkAEfkw8E/AL+FkJitIgANXAcdU9dnodfcDbwASBfgbgC+rE23msIh0iMh5qvpCgffMiVIIgnzMWfW+3hWb8CycXSASdEJURkIRZE4IB8KcefTMqnPKMQmp93Y0iiNZ6DV5m+KfSyXc8l3CSlQGEjOi+fy+eJ6CfCfqqdbhTz10ir4dfTRwbmJhfX/tkzUSm4g8AVyhqoHo5wbgEVW9VER+oqovK+jGIm8BXq2q74h+vgV4haruTSjzTeBjqvr96OcHgD9R1VVh1kTkNuA2gL6+viv3799fSLUAGH9iHJfXWV1wdbiITDnO95FghJ5Le7KeH5gLMHNyBnEJ4pZ4Tuy2LW2rzG8Ak89OEglGEHdCBKiw4vK66Dy/NvYbZiJW/1QZlxo3NxIYCdC3o69i9ajXdkzF3NwcLS0t2Qsaq4i1XeL7nEiu73MygbkA82PzhJfDuBvc+Hv9Kd9rWNknA3MBRKJ9U5wEK4XUI1U/D84HAfD6vfFjxfZ963uFU+q2u+6661JGYstFA/97HO3369HPNwL/ICJ+VmrL+ZIqKEzybCKXMs5B1X3APnBCqRYTxu7Apw8wMzZDQ1sDDW9oYPnry/HweYO/l/266UIiaq+mDL83LOdm1KFAiIljEwQXgmzauYmX3fGymjeB7XufE4Jy+NAwoaWQE4oymov84v99MeOfGeetT7617PVIbMdShE+tBdZzOMtiibVd4vscI5/3GQp3EIu9G+ISTj98mvByGJfHRWg5xLZXbsu7HsnXjDF/xonZvWnnppL1fet7hVOptsvqxKaqfwq8E5gCpoF3qeodqjqvqjcXce/ngS0JnzcDpwsoU3KKzSqTztFl5NGRlI40MdMfAqNHRxGEvpc4Gms9OKLEnIca2hrwNnsRl8Rn/26fuyLaN9S2g59RPYp1mhw+NMx977yP4aFhJp6eYHgo+jnhvUznJJfoWNexrYNIOEJoKYSnycPyzDKzI7PMj82ndGo7fOdhPrn5k/yZ/8/45OZPcvjOw6uuGcPT4GHTzk3W99cZOW0jU9WHcda7S8kR4CIR2Q6cAm4CfiOpzDeAvdH18VcA0+Ve/4bs+5gh84w8laPL9PPTLJxZWOUlmphBzL/Bz6adm1acB7XviBJbF2ze0ExwPoinwYM2KB0DHYgIO/eusvyUjVp08DOqS7FOk0MfGmL2hVk8Pg/uBjcaVmZfmGXoQ0PceujWjHvDE9fMm7qaaN/WzvSJaRrbGh37YtSemHzeyCMjPPDBB3B73PH18wc++ACQfh3eBPb6o2rpRFU1BOwFvgU8Afyjqj4mIu8SkXdFi90PPAscAz4HvLtS9Ytt/ei5tGdVRKZsWzZSzfhjCQ0yRVaq14QZifHJ/X1+fK0+Wvtb6bm4h7YtbTaoGFWnmMQeI0dH8Pg8uDwuRASXx4XH52Hk6AiQOWpaslWo5+IebvraTex9ci/+DX5az2tNed6DH38Qt8eNp9GDiOBp9OD2uHnw4w+apcmIU9VALqp6P46QTjz22YT/K/D7la5XNrJ5taaa8ft7/bRvbl9xnWThXM4AJ+UmneY7NDRU8boYRikRhGRnX1VFoi462XY/pHs3Mp23OLno5DRPwN3gZnFyMeM1jfWFRWIrgFy2KyW/YKkc25KFczmirBmGURx9V/Rx+uHTce07EooQDobZeOVGoPCJd6bz5kbmCMwG8DSeG6ID8wGIOPnPY8t2QEmir6VaEjRqn6qZ0OuZQhIU5OJIY6ax+qCU+diN2mfwI4O0bmwFFwSXguCC1o2tDH5kECjcSS7Tede8/xrCoTChpRCqyvLcMqHFEP6N/viy3X233cd977yv6Ohr6ZYEk8c4o/YwAV4AhbywuQrnQsIuGpXDQlauPwZ2D3DjvhvZPrid7ou62T64nRv33Rh/dwudeGc67+r3XM31H70eX6svrnm3b2+n77K++Hr54sQic6NznH36LCd+cIKzT58lHAznnbEs3Rr+/Nh8gS1mVAozoRdAoV6t+axb5Zv1qNqsl0QKFq51fZLt3S10TTrTeVe/5+p41sJ9O/etWrYLzAeILEdwe9y4fW7Cy2Gmh6cJLYbyqkO6JcHwcjiv6xiVxwR4gZTbiaSeBEWmycZaw8K1GtUg1Xq5hhVc4PI4hlTxCJFQhOBCsOhrB+YCNDU0labyRtkwE3qNUuyWskzrtKVew820jWatUYj/g2EUS6plOxHB7XMTCUVDPYciKLrKe72QawfmA/h7/eV4FKOEmACvUYoRFJnWacuxhluv+9cLoVypUA0jE6nWyzdfvZmuC7twN7gJLYdwN7jpGOjIO/JhurX4dPHdjdrBTOg1SrYtZZnWnDOZ34GMpvlC1rLref96vpQrFaphZCN52S42Ge++qHvFGFHIZDLVkuDw0HBxFTbKjmngNUomD9VsWnQmjTjTd4Vq5+tNKy0mqpdhFEry0hdg207XOaaB1zDpHOWyObhl04jTfZfuukMfHsK/wZ9WK8+klSbO4teLp7phlJpMjqKpMhwa6wMT4HVIzBN6YXyBqeNTBBeCeJu8zL4wC2Q3v6f77uB7D67ysA4FQoweddIUZtrOls0rv962xRm1y3qcCNbTrhSjcpgJvcoU4hHu8/s48YMTnPrRKRbGF1CU4EKQhTMLDB8azmh+z/RdKse5iWMT+JqL9zBfT57qRvlYr4F01pOjqJE7poFXkUK00uFDw8yOzDqCNpqOMDgXdDIkNXjY/8b9XPCqC9i1d1da01qis9vkc5NxIZpKcw8uBOO5yWMUMnDY/mmjFKxXTXQ9OYoauWMaeBUpRCs9ctcRWvtb8TR64gEcEGcPqLgFImTVStJpMbDaKWbTzk14GlbO8woZOKq1f9rilq8t1pommmv/XG+OokZumACvIoUMRrFzYsK+qasJcTlpDUUEr9+bdSKQLX9xoof14B2DJRk4qjEArVdz61qmngLpZBPO+fRPS3RkpMIEeBVJHowWxhc49aNTTD4zmXY2HjunY6CDSDhCJBQhEo4gLiESjtCxrQPIPBHIZ+JQqoGjGgOQrbuvPepFE81FOOfbP237opGMrYFXkcQ159ByiLHHxgDovbw37Xp47Byf30fPpT1MHJtAEDxNHjZcuoHmnmYgs1aS73paqeK+lzt+fDK27l57FOtBXi+BdHJZq7f+aRSLCfAqkjgYPfOdZ/A0eui+qJvm7uZ4mWTnnOQBbPvgdrZcu4Wjdx/F7XOjEV21bSyZbNvM1grm+FMbxIT2yKMjLJxZoG1rG+2b2wveSljpiWAh5CKcK9U/1+O2u/WCCfAqExuMYukCY+vZkNmsnfwC9r+0P2etpF60mHTkOiCtl4lKLZO40yIwG0AjyvTxaXx+X3yiuhY9yHMRzpXonxZ/YW1jArxGiL3w4WCYqWEnOIvL66L3xb05nZ+vVlIPWkwq8hmQkicqvhYfPr+Pg+89aJpIhUg0JQcXg7gb3GhYmRqeorm7ec2ajHMRzpWYSK/XbXfrBXNiqxF27d3F7Mgsoz8dJbQUwuV2EVoKMXt61rymEyjU8eeGT9wQdxg0j/TKkegw6W32omHF5XHFc1av1SWNXJ02y+2Ytta23RkrMQ28RhjYPUBrfytLE0uEg2E8zR66L+7G7XPbbDmBQh1/TBOpDomm5I5tHZx54gyRUARPsyfuQb5WlzRqwcplfiBrGxPgNURgPsCmqzatWAfXiNpsOYFUA9L089Msnl1k3859aU3j5vFbHRJNyU1dTbRva2f6xDSNbY34e/1V870oxLGrGs5gxd7T/EDWNlUxoYtIl4h8W0Sejv67ajooIltE5D9E5AkReUxE3lONulaSegpSkQ+ljIaWvA946sQUE8cmaOpuymgaX6ttW+skm5J7Lu7hpq/dxN4n91ZtL3MhAX6qERSoFPdMbn9wBPnB9x60yIRrgGqtgX8AeEBVLwIeiH5OJgS8V1UvBa4Gfl9ELqtgHStOvQSpyIfAXKCkA1/ygLR4dpGuC7vo2NqRcU18LbZtvVCtACSBuUDKiWOhIYwrHRSoVPdc4QcyHwA1P5C1QrUE+BuAL0X//yXgjckFVPUFVf1x9P+zwBPApkpVsBqsxXCJ82PzJR/4EgVCS38L7ZvbV3yfyjS+FtvWSM/woWFmTs6knDgWE8I4n3OKpdT3tMiEaw9R1crfVGRKVTsSPk+qalpbpogMAN8FLlfVmTRlbgNuA+jr67ty//79Janr3NwcLS0t2QsaKZk8M4lOr+5jkWCEnkt7ir/+s5NEgtFELlE0rLi8LjrPr3/zuPW/7ATmAs4WzOUw7gY3/l4/82PzSJsQmYrEy8X6BZB3n6lGPyv1PcefGI8/fyLp3kXre4VT6ra77rrrHlbVncnHy+bEJiLfAfpTfPXf87xOC/BV4P9NJ7wBVHUfsA9g586dOjg4mM9t0jI0NESprrUe+ecv/DMzd8+scDpbnll2HJh+b7Do6w/LuX3hiU46a0W7tv53jlQOXQAHb1/5+5+ZP8Py7DKb3r+J5a8vx8/XiDI3OscNn7gh7z5TjX5W6nse+PQBZsZyfxet7xVOpdqubCZ0VX2Vql6e4u/rwKiInAcQ/Xcs1TVExIsjvL+iqveWq65G+fD3+su69mym8fVBOoeuoQ8PpTQLB+eDaHil5SfmtFhInym2nxXiyFnqvm1+IGuPam0j+wbwW8DHov9+PbmAiAjwN8ATqvrJylbPSKSYrSy+Fl/KaFMAB/YcKMmWnFrYb2uUl3T7+EcfHWXLNVtWlPW1+JygMVEhlWr7VCF9ptB+Vkw401L27XoPoWysploC/GPAP4rI7wAngD0AIrIR+Lyqvhb4ReAW4Kci8kj0vA+q6v1VqO+6pRSxlJMHIYvPbORLun38ipO8JzlQSf8V/bRtaUN7terCqpaCCNlkd21RFQGuqmeB61McPw28Nvr/7wOSXMaoLOUYfGppQDNqh0yWnnQRxfp39Dtbo1gdqGRYh9lzYE81HmUFFkTIKBcWC93ISDm2zxRzzVIGhTFqh2xBS9Kt3w7eMVjzPhAWRMgoFxZK1chIOWIpF3pNM72vXbJZZbKt39ZySFMLZ2qUC9PAjYyUw3O10GtaIIq1Sy5WmWpFdCs2pKntlDDKhWngRkbK4bla6DVtLXHtUstZs0rhs2HOY0Y5MAG+DsnXHFiOwSfTNdPVr5YHeaM4atnMXMqJYzUymhlrFzOhrzOqkVUpHzLVzwJRrF1q2cxcKie0Wn/3jPrDNPAaphyz9VrfwpWpfrGMShaIYm1Sq2bmUlkHav3dM+oPE+A1Srk8rmt9HTlb/Wp1kDfWLqXyA6n1d8+oP0yA1yjlmq3X+jpyrdfPWJ+UYuJofdsoNbYGXqOUK/9wra8j13r9DKNQrG8bpcYEeI1SruhNtewsBLVfP8MoFOvbRqkxE3qNUs5tNbW+jlzr9TOMQrG+bZQSE+A1ylpP/Wf7YQ3DMIrDBHgNs1Zn6xbT3DAMo3hMgBslJ1G73vDuDQzL8ArBbPthjbWMWZeMSmFObEZJSY42FQlGVkWbSuVhH1oO8cx3nrE0oUZdY9HWjEpiAtwoKTHtOhwI88JPXiC4GGT6+DRDHxqKl0n2sF8YX2DssTFcHpcNekZdYxnzjEpiAtwoKZPPTRJaDnHmiTOEl8OICBpRTj10Ki6Qk/fDThybAKD7om4b9Iy6ppD4DcOHhjmw54BZn4y8MQFulJTO7Z1MPjOJy+3C5TnXvbzN3rhATt4PGwlF6L28l+bu5nh5CzFp1CP5xm8wk7tRDCbAjZKya+8uAgsBVDV+LBKO0HVh1wqBPLB7gD0H9nDbQ7dx/qvOx+Nb6U9pISaNWiad1pxvtDUzuRvFYALcKCkDuwfYtGsT4hZCyyEQ2HDZBjw+T1qBbCEmjXoik9acb7S1coVMNtYHto3MKDmDHxmM7/P2Nntxe90Zo8it9aA1xtoi2zbIfOI3WIIToxhMgBslJ1EgR4IR/L3+rAJ5rQatMdYepUwLWs6QycbapyomdBHpEpFvi8jT0X/TTjdFxC0iPxGRb1ayjkZxxNa4ey7tYc+BPSacjTVDKRMNWYIToxiqpYF/AHhAVT8mIh+Ifv6TNGXfAzwBtFWqcrWERXUyjNqi1FqzWZ+MQqmWE9sbgC9F//8l4I2pConIZuB1wOcrU63awraYGEbtYVqzUStUSwPvU9UXAFT1BRHpTVPuU8AfA62VqlgtYTHDDaM2Ma3ZqAUkcb9uSS8s8h2gP8VX/x34kqp2JJSdVNUVC0gi8qvAa1X13SIyCLxPVX81w/1uA24D6Ovru3L//v1FPwPA3NwcLS0t2QuWgfEnxnF5VxtJIsEIPZf2VKFGDoG5APNj84SXw7gb3Ph7/au2wsSItV8+5xjnqGb/q3es7YrD2q9wSt1211133cOqujP5eNkEeCZE5ElgMKp9nwcMqeqLksr8L+AWIAQ04qyB36uqv5nt+jt37tSHHnqoJHUdGhpicHCwJNfKlwN7DqzaYrI8s4y/18+eA3uqUqfEVKCJ63/pTIhDQ0MMyEBe5xjnqGb/q3es7YrD2q9wSt12IpJSgFdrDfwbwG9F//9bwNeTC6jqf1PVzao6ANwE/HsuwnstUYsBTgqJHGXRpgzDMEpPtQT4x4BfEZGngV+JfkZENorI/VWqU81Ri84yhUSOsmhThmEYpacqTmyqeha4PsXx08BrUxwfAobKXrEapNacZRIjRy2ML3D2qbMsTS8hLuGLu7/I4B2rA7ZYtCnDMIzSY7HQjbyImfWnTkwxenSUxalFVBVPo4fTPz7Nfe+8b9U2t1pcCjAMw6h3TIBXiXrNARwz6y+eXSS0HMLldtHY1oi3yYvH52FpcmnV2nYtLgUYhmHUOxYLvQokenInBmgpVqhVImpb7B6Lk4sA+Pw+3D43AOIWwsvhlGvbtbYUYBiGUe+YBl4FyuGVXcqobemsA4n3aGxrBIHl2WXCwTAAGlbcPretbRuGYVQAE+BVoBxe2aWaFGSaCCTeo2OgA4/Pg6qyPLtMJBQhFAjR2Nloa9uGYRgVwAR4FShlNqMYpZoUZJoIJN6juaeZvh19NLY3gkIkEmHjyzdy4+duNFO5YRhGBbA18CpQjhzApdqqlSnXcfI9mnua6fP1VTUynGEYxnrFNPAqUA6v7FJt1cpkHbDtYIZhGLWDaeBVotRe2bFJQaIX+uDe1UFVspHJOlCqexiGYRjFYwJ8DVGKSUE2IW3bwQzDMGoDE+DGKkxIG4Zh1D62Bm4YhmEYdYgJcMMwDMOoQ0yAG4ZhGEYdYgLcMAzDMOoQE+CGYRiGUYeIqla7DiVHRM4Ax0t0uR5gvETXWo9Y+xWHtV/hWNsVh7Vf4ZS67bap6obkg2tSgJcSEXlIVXdWux71irVfcVj7FY61XXFY+xVOpdrOTOiGYRiGUYeYADcMwzCMOsQEeHb2VbsCdY61X3FY+xWOtV1xWPsVTkXaztbADcMwDKMOMQ3cMAzDMOoQE+CGYRiGUYeYAI8iIq8WkSdF5JiIfCDF9yIifxn9/qiIvLwa9axVcmi/m6PtdlREHhSRK6pRz1okW9sllNslImEReUsl61fr5NJ+IjIoIo+IyGMicqjSdaxVcnhv20XkPhF5NNp2b69GPWsVEfmCiIyJyM/SfF9euaGq6/4PcAPPAOcDPuBR4LKkMq8F/hUQ4GrgP6td71r5y7H9rgE6o/9/jbVf7m2XUO7fgfuBt1S73rXyl2Pf6wAeB7ZGP/dWu9618Jdj230Q+PPo/zcAE4Cv2nWvlT/gl4CXAz9L831Z5YZp4A5XAcdU9VlVDQD7gTcklXkD8GV1OAx0iMh5la5ojZK1/VT1QVWdjH48DGyucB1rlVz6HsAfAF8FxipZuTogl/b7DeBeVT0BoKrWhg65tJ0CrSIiQAuOAA9Vtpq1i6p+F6dN0lFWuWEC3GETcDLh8/PRY/mWWa/k2za/gzMrNXJoOxHZBLwJ+GwF61Uv5NL3LgY6RWRIRB4WkbdVrHa1TS5tdxdwKXAa+CnwHlWNVKZ6a4Kyyg1PqS5U50iKY8n763Ips17JuW1E5DocAX5tWWtUP+TSdp8C/kRVw44iZCSQS/t5gCuB64Em4IciclhVnyp35WqcXNruvwCPAL8MXAB8W0S+p6ozZa7bWqGscsMEuMPzwJaEz5txZpz5llmv5NQ2IrID+DzwGlU9W6G61Tq5tN1OYH9UePcArxWRkKp+rSI1rG1yfXfHVXUemBeR7wJXAOtdgOfSdm8HPqbOgu4xEXkOuAT4UWWqWPeUVW6YCd3hCHCRiGwXER9wE/CNpDLfAN4W9Sq8GphW1RcqXdEaJWv7ichW4F7gFtN8VpC17VR1u6oOqOoA8E/Au014x8nl3f068EoR8YhIM/AK4IkK17MWyaXtTuBYLhCRPuBFwLMVrWV9U1a5YRo4oKohEdkLfAvHM/MLqvqYiLwr+v1ncbx/XwscAxZwZqYGObffh4Bu4K+jmmRILdNRrm1npCGX9lPVJ0Tk34CjQAT4vKqm3Paznsix7/0p8EUR+SmOOfhPVNVSjEYRkX8ABoEeEXke+DDghcrIDQulahiGYRh1iJnQDcMwDKMOMQFuGIZhGHWICXDDMAzDqENMgBuGYRhGHWIC3DAMwzDqEBPghpGGaOavR0TkZyJyQESaRWQgXeahMtWhSUQOiYg7+vnj0axQH89wzu0i8r4Ux/Ouu4h8MVX2s2hY0oK2AeaaoUlE/iaaBeuoiPyTiLREj3eKyD9Hj/9IRC5POs8tIj8RkW/mWJ+5Qp4jw/W+IyKdpbymYaTCBLhhpGdRVV+qqpcDAeBdVajDb+Mk4ghHP/8u8HJVfX8V6lIqXgNcFP27DfhMmnL/VVWvUNUdOAFF9kaPfxB4JHr8bcCdSee9h+oGarkbeHcV72+sE0yAG0ZufA+4MPp/t4h8LqoJHxSRJgAReaeIHIlqjV+NRv1CRPZEtfhHo2E8Y1rix6Plj4rI76a57804kcQQkW8AfuA/ReStIrJNRB6Inv9ANNrdCkTkyuh9fwj8fsLxlPePasd3icjjIvIvQG+GNvlNcXK7/0xErhIRl4g8LSIbotdyRbXsnqTzcsrQFIu3LU7knybOxZC+DHggWubnwEA0Shgishl4HU7I3pSIE3nsh9Fn/9OE4y3RdvyxiPxURN4QPf6nIvKehHJ/JiJ/KCLnich3E6w0r4wW+Qbw6xnazTBKgglww8iCiHhwtMafRg9dBHxaVV8MTAG/Fj1+r6ruUtUrcDTA34ke/xDwX6LHXx899js4YRV3AbuAd4rI9qT7+oDzVXUYQFVfzzmrwD04maK+HNVEvwL8ZYrq/y3wh6r6C0nH093/TTjhMl8CvBMnj3s6/Kp6DY62+YVolqq/w5l0ALwKeDRF5K6cMzSJyN8CIzjxt/8qevhR4M3R768CtnEuPe2ngD/GibiWjjuBz0SffSTh+BLwJlV9OXAd8Ino5OFvgN+K3s+FE3L0KzhpSr+lqi/Fia3+CEA0bW6DiHRnqINhFI0JcMNIT5OIPAI8hGPC/Zvo8edU9ZHo/x8GBqL/v1xEvidO2MmbgRdHj/8AJxzlO3FCVgLcgBMj+RHgP3HCzF6UdP8enAlCOn4B+Pvo/+8mKcObiLQDHap6KKFMjHT3/yXgH1Q1rKqngX/PcP9/gHhO5DYR6QC+gGPWBsf8/7cpzss5Q5Oqvh3YiDMhemv08Mdw0oM+gpMn/SdASER+FRhT1Ycz1BngF2N1Z2WbCPBRETkKfAdnUtEXnUCdFZGX4bTbT6LJeI4AbxeR24GXqOpswrXGovU2jLJhsdANIz2LUe0qjqOQsZxwKIxj3gX4IvBGVX1URG7FiZGMqr5LRF6BY9p9REReiiMs/kBVv5Xp/kBjHvVNlQI3XazklPcXkddmOCfb/VRVT4rIqIj8Mk7SkJtTnJdXhqZoGtV7gPcDfxs1rb89Wl8Bnov+3QS8PvoMjTiTir9T1d/Moe5E67oBuFJVgyIyzLn2/zxwK9CPM0lBVb8rIr+E87veLSIfV9UvR8s34vx+hlE2TAM3jNLRCrwgIl4SBJeIXKCq/6mqHwLGcYTXt4Dfi5ZFRC4WEX/ixaKmWLeIpBPiD+IILaL3+37S+VPAtIhcm1AmRrr7fxe4KbpGfh6OKTkdb42eey2OOX46evzzOKb0f0xwvkska4am6HcXxv4P3Aj8PPq5I7q8APAO4LuqOqOq/01VN0eztt0E/Hsa4f0DVrZbjHYcDT4oTt76bQnf/TPwapzlhm9F67EtWv5zONaZlyfUtx8YTtVohlEqTAM3jNLxP3HM0cdx1stbo8c/LiIX4Wi9D+Cs4R7FMb3/ODrgnwHemOKaB3FM499J8d0fAl8QkfdHz0+V6ejt0TILRAVPlM+nuf8/A78crf9TwCHSMykiDwJtOObyGN/AMZ2nMp9DhgxNInI/jlAeAb4kIm047fYo8HvRYpcCXxaRMPA453wNcuU9wN9HHdO+mnD8K8B9IvIQznr2z2NfqGpARP4DmEqYlAwC7xeRIDDHuaWDK4HDqhrKs16GkReWjcwwapjouusfqeot1a5LroizP/z/qOorsxauE6LOaz8G9qjq01nK3gl8Q1UfqEjljHWLmdANo4ZR1Z8A/yHRQC61joh8AEer/W/VrkupEJHLcKwFD2QT3lF+ZsLbqASmgRuGYRhGHWIauGEYhmHUISbADcMwDKMOMQFuGIZhGHWICXDDMAzDqENMgBuGYRhGHfL/A38fXDoHJ6FWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 确保两波段在相同时间附近有观测，按最近配对 g 和 r 波段的观测点\n",
    "from scipy.spatial import cKDTree\n",
    "\n",
    "# 构建 MJD 和 calmag 映射（仅保留必要列）\n",
    "r_tbl = df_r[['mjd', 'calmag']].dropna().reset_index(drop=True)\n",
    "g_tbl = df_g[['mjd', 'calmag']].dropna().reset_index(drop=True)\n",
    "\n",
    "# 用 KDTree 匹配时间最接近的点\n",
    "tree_r = cKDTree(r_tbl[['mjd']])\n",
    "dists, idxs = tree_r.query(g_tbl[['mjd']], distance_upper_bound=0.01)  # 匹配时间差在0.01天内\n",
    "\n",
    "# 筛选有效匹配点\n",
    "valid = idxs < len(r_tbl)\n",
    "matched_g = g_tbl[valid].copy()\n",
    "matched_r = r_tbl.iloc[idxs[valid]].reset_index(drop=True)\n",
    "\n",
    "# 计算 g - r 颜色\n",
    "matched_g['rmag'] = matched_r['calmag']\n",
    "matched_g['phase'] = (matched_g['mjd'] / true_period) % 1\n",
    "matched_g['g-r'] = matched_g['calmag'] - matched_g['rmag']\n",
    "\n",
    "# 按相位绘图\n",
    "plt.figure(figsize=(7, 4))\n",
    "plt.scatter(matched_g['phase'], matched_g['g-r'], color='purple', alpha=0.7)\n",
    "plt.xlabel(\"Phase (folded by 0.394 days)\")\n",
    "plt.ylabel(\"g - r Color\")\n",
    "plt.title(\"Color Variation with Phase\")\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6048a05a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
